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Abstract. We calculate analytic expressions for the distribution of bipartite 
entanglement for pure random quantum states. All moments of the purity 
distribution are derived and an asymptotic expansion for the distribution itself 
is deduced. An approximate expression for moments and distribution of Meyer- 
Wallach entanglement for random pure states is then obtained. 
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Introduction 

The question of generating and measuring entanglement in multipartite quantum 
systems has become of greater interest with the development of the field of quantum 
information. Entanglement generation is an important aspect of several quantum 
information processes, such as superdense coding [1], quantum communication [2], or 
quantum data hiding [3]. Various methods have been proposed in order to generate 
highly entangled quantum states, based on pseudo-random unitary operators [4] or on 
the entangling power of chaotic quantum maps [5, 6] or intermediate quantum maps 
[7] . Entanglement generation by means of pseudo-random unitary operators or chaotic 
quantum maps relies upon the fact that unitary evolution of any initial state leads to 
states whose entanglement properties are close to those of random states, in particular 
to highly entangled states. 

In order to quantify the entanglement of a state, or the entangling power of an 
operator, a number of entanglement measures have been proposed, based either on 
quantum information theory or on thermodynamical considerations: entanglement of 
formation and distillable entanglement [8], relative entropy [9, 10, 11], n-tangle [12], 
concurrence [13]. For bipartite entanglement of pure states, these measures all reduce 
to the entropy of entanglement [14], which can be proved to be a unique entanglement 
measure [15, 16]. The entropy of entanglement corresponds to the von Neumann 
entropy of the partial density matrix obtained by tracing over one subsystem. Rather 
than the von Neumann entropy itself, one often prefers to consider the purity R, which 
corresponds (up to constants) to the so-called linear entropy, that is the first-order 
term in the expansion of the von Neumann entropy around its maximum. To quantify 
the degree of entanglement of multipartite pure states, one measure commonly used, 
based on purity, is the measure proposed by Meyer and Wallach in [17]. It consists in 
taking the average of the bipartite entanglement of one qubit with all others, measured 
by the purity (see Equation (18)) [18]. Meyer- Wallach entanglement was used e.g. to 
quantify entanglement generation for pseudo-random operators [4] or intermediate or 
chaotic quantum maps [7, 19, 20]. 

The study of purity or Meyer- Wallach entanglement is of particular interest for random 
quantum states. Random pure states as column vectors of random unitary matrices 
distributed according to the invariant Haar measure can be shown to be entangled with 
high probability. Various analytical calculations have been carried out to characterize 
entanglement properties of random states. Expressions for the first moment of the 
purity have been obtained by Lubkin [21]; the second and third moments have been 
derived in [6], following earlier work [23]. The average entropy has been obtained 
in [22]. Statistical properties of entanglement measures for random density matrices 
were obtained in [25]- [27]. The average value for each Schmidt coefficient of a random 
pure state has been calculated in [28] . In [29] , the average entropy of a subsystem was 
obtained from the average Tsallis entropy [30]. 

To further characterize entanglement of random pure states, our aim here is to give 
an exact expression for all moments of the probability density distribution P{R) of 
the purity for a bipartite random pure state. Since the probability distribution P{R) 
is defined over a bounded interval {R is bounded), the knowledge of all moments 
determines uniquely the probability distribution [31]. There are various techniques 
to obtain a function approximating the exact probability density distribution in a 
controlled way (that is, by an expansion where the error can be bounded) from the 
knowledge of its moments. In [32] an algorithm was given to construct polynomials 
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converging to the probability distribution. We will rather follow [33], where the 
asymptotic expansion for nearly gaussian distributions is calculated at all orders. 
In Section 1 the moments (i?") for the distribution P{R) are calculated, and the 
construction of the asymptotic expansion of P{R) at all orders from its moments 
is recalled. In Section 2 the approximate moments (Q") for the distribution P{Q) 
are derived. For both distributions, the moments are expressed as sums involving a 
finite number of combinatorial terms and can be easily calculated effectively. As an 
illustration, we give the first values of the cumulants and calculate the probability 
density distribution expansion for Meyer- Wallach entanglement. 



1. Bipartite entanglement for random pure states 

Let ^ he a pure state belonging to a Hilbert space Ha Wb, where Ha and Hb are 
spanned respectively by {\ai)}i<i<p and {|6;)}i<i<q. We assume that p < q. Let Xi 
be the Schmidt coefficients for "if. That is, the state "if has a Schmidt decomposition 
(sec e.g. [34]) 

p 

I*) = X)v^l«i>® l^i)- (1) 

i=l 

The bipartite entanglement measure for ^ can be expressed throught Schmidt 

coefficients Xi. The entropy of entanglement is the Shannon entropy of the Xi's: 
S{'^) = — J2i=i lna;i. The purity i?(\I') of the state \1/ can be expressed as 

i?(*) = £.r2. (2) 

i=l 

For random states the Schmidt coefficients are distributed according to the density 

p{xu...,xp)=j^ n i^i-^jf n 4"'^fi-E^^) (3) 

l<i<j<p i<k<p \ i=l / 

for Xi e [0, 1], with some normalisation factor Af [37, 6]. The n-th moment of the 
purity is then given by 



(ii") =Af [ dxi...dxp Y[ {xi - Xjf Y[ xl 



Q-p 

l<i<j<p l<k<p 

■P 



(4) 



(5) 



X {xl+xl + --- + '^^ (^1 - E j • 

The calculation of (i?") requires the evaluation of integrals of the form 

7(n) = J dxi... dXpV{-x.fx{ . . . a;^^ ^1 - ^ Xi^ /n(x), 

where r = q — p, x = {xi, . . . ,Xp) and n = (ni,...,np). The function /n(x) = 
{x"^X2^ . . .Xp" + all permutations of the rij} is a symmetric function of the Xi, and 
V is the Vandermonde determinant 

F(x)= n i^i-^i)- (6) 

1<2<J<P 
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The integral I{n) is evaluated in the appendix and yields 

where " +pcrm" indicates that the expression (7) is a sum over all permutations of the 
Hi. Now the function / for a given (ii") is obtained by multinomial expansion of the 
term 

ni+n2-\ \-np=n ^ 

The normalization constant J\f in (4) is given by the choice n = in (5), i.e. 
^= l//(0). This leads to 

/^nv ^ (p'+rp-l)! ^ n\ 

^ ' (»2 +rp + 2n- 1)! ^ ni!n2!...n„! 

nl+r^2^ Vnp=n ' 

., nLi('^ + 2ni + z-l)! T-r 2n,--2n»+j-i 

Replacing r by its value q — p and correspondingly changing all indices i to p + 1 — i 
(and rii to rip+i-i), one finally obtains 



(pg + 2n - 1)! nilns! . . . n^! 



><nn^^W^ n (2n.-.-2n,+,). (10) 
i=i W i<i<j<p 

Note that one can cast (10) into an expression more symmetric in p and q by noting 
that 



ni<j(2»»--<-2nj+j) 



n 



nLi(?> + 2n,-i)! 



2n,- 



(2n,)!Ai^ 2ni+i 



,(11) 



yielding 



im\ = (Pg- 1)' 

^ ' (pg + 2n-l)! ^ f^^ mlns! . . . rip! ^ ^ 



xn 



{q + 2nj - i)!(p + 2nj - i)! /^i _ ^"-J' 



(g-i)!(p-i)!(2ni)! V 2ni+i-i 



Equation (10) is a closed expression, involving only a finite sum over partitions of n 
into numbers greater or equal to 0. Note that the order of the ni matters: for instance 
for p = 2 and n = 2 the sum will involve three terms (ni, 712) = (2, 0), (1, 1) and (0, 2). 
These partitions can be easily generated for any n by some suitable algorithm (see 
e.g. [33] for such an algorithm generating the partitions required). From Equation 
(10) one can get the expressions for the cumulants of the distribution P{R}- Indeed, 
given the moments /i„ of a distribution the n-th cumulant reads (see e.g. [33]) 

»» = "'E(-ir-('-i)!ni(S)"'- w 

{fe™} rn=l "* 
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where r = fci + • • • + fc„, and the sum runs over all fcj > 0, 1 < i < n such that 
kl + 2^2 + ... + nkn = n. As an example, the first five cumulants read 

M = (14) 
1+pq 



K2 = 
K3 = 
K4 

K5 



2(p2 _ l)(g2 _ 1) 

(1+P(Z)2(2+m)(3+m) 

8(p2-l)(g2-l)(p + g)(-5+pg) 
(1 + pg)3 (2 + pq) (3 + pq) (4 + m) (5 + pq) 

{l+pq)^{2+pq){3+pq) ULiii + PQ) 

(1 + pq)^{2 + pq){3 + pq) ULi + Pi) 



where Ap g and Bp q are polynomials in p and g defined by Ap q = 28 — 112^^ — 153pg — 
79p^q-li2q^ -98p^q^ -Up'^q^ -79pq^ -3p^q^ +p^q^ -llp^q*+Ap*q^+p^q^ and Bp,, = 
3528 - 6552p2 _ 6343pg - U9p^q - 6552g2 + 1545p2g2 + i237/g2 _ 449^^3 _^ 1164^3^3 + 
132p5r73^;L237p2g4_274pV-41pV + 132p35''^-93F'^g''^+p^g^-41p4g6_^9/g6_|.p5^7 

As expected, ki corresponds to Lubkin's expression [21] for the average purity. For 
n = 2, 3 one recovers the expressions derived in [6]. For larger n it is easy to generate 
the exact value for each cumulant. 

In the case p = 2, the analytic expression for the probability distribution P{R) dR 
can easily be obtained analytically directly from (2)- (3). It reads 

P{R) dR = A{1 - i?)«-V2ii- 1 dR (15) 

for 1/2 < R < 1, otherwise {A is the normalization factor). For p > 3, the 
asymptotic expansion of the distribution can be obtained (see [33] and references 
therein) by Edgeworth expansion as a function of the normal distribution Z{x) = 
exp(— a;2/2)/\/27r, the mean = ki, the variance cr^ = — /if = K2, and rescaled 
cumulants 7r = Hr/cr^^^^: 

P{R) = \Z [1+ (16) 



u 



V - / ,km\ \{m + 2)\ 

" — ^ {km} ^ "1=1 ' 

For each s the sum runs over kj > such that ^jjkj = s, and t is defined by 
t = Y^jkj. The IIe„(a;) are Chebyshev-Hermite polynomials defined by IIe„(a;) = 
(— l)"e^^/25n(;-2;^/2 (^iierc d is the differential operator with respect to x) and 
correspond to rescaled Hermitc polynomials: 

["/^l (_'\\k^n-2k 

fe=o ^ ' 

Equations (13)-(17) together with the knowledge of the moments (10) allow to obtain 
explicitely the asymptotic expansion of the probability density distribution at any 
order. 
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2. Multipartite entanglement 

The Meyer- Wallach entanglement of a pure M-dimensional state 4* coded on m qubits 
(with M = 2™) can be defined by 

QW=2(^l-^E^fc)> (18) 

where Rk is the purity (2) of the fc-th qubit [18]. In order to calculate Q{'^) for 
bipartite random pure states we need to calculate the average purity of a bipartite 
system belonging to a Hilbert space Ha (^Hb, where Ha has dimension p = 2 and 
Hb bas dimension q = 2"^~^ = M/2. The moments (i?") can be obtained in this case 
either from Equation (10) or directly from the distribution (15). In both cases it leads 
to 

r(3+i)^/n\ (fc + i)! 



The calculation of the moments (Q") involves terms of the form ((X^j i?i)*')- These 
terms depend on correlations between the purities i?;. However if we make the 
assumption that for two different qubits i ^ j we have (RiRj) = {Ri){Rj), we get 
a distribution P{Q) which turns out to be very close to the numerical distribution 
obtained by generating random matrices. Making this assumption we get 



The n-th moment is then 



Gathering together terms having the same exponents, we finally get 

(Q") = 2" ^ (^^^ ^-1— ^ ,,,rk\{m~ r)\ H (^) '^^^^ 

fe=0 ^ ' {rk} *=1 

where r = ^^'^ (-^") given by (19). The sum runs over all > such that 
^jjrj = k. Prom Equations (19) and (22) one can now obtain the cumulants for the 
distribution P{Q) for an m-qubit system (M = 2™). The first ones read 

= ^ (23) 
6(M-2) 
2 (M + l)2(M + 3)m 

24(-M2 + 7M - 10) 



/to — 



(M + 1)3(M + 3)(M + 5)m2 
Q _ 144(M4 - 12M3 + QNP + 133M - 210) 
^ ~ (M + 1)4(M + 3)2(M + 5)(M + 7)m3 
Q _ 1152(1890 - 1763M + 337M2 + 78M3 ~ 23M4 + M^) 
^ ~ (M + 1)5(M + 3)2(M + 5)(M + 7)(M + %)m^ " 
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We can use these approximate cumulants to obtain an analytical formula for P{Q). 
Calculating the first terms in the asymptotic expansion (16) we obtain (the first terms 
can be found in [38]) 

Q - 



P(Q)^iz(^ 



1 + ^Hes 
6 



_He4(^ 



+ 



T5„ (Q- , '^3T4„ 



a 



1296 



Heg 



(24) 
(25) 



wit 1 n = nf, 



and Ti = Ki/a'^, i > 3. Figure 1 displays the probability 

density function P{Q) for m = 10 qubits as obtained by averaging over numerically 
generated random matrices, together with the plot of analytical expression (24) trun- 
cated at order (gaussian), 1 (first line of (24)), 2 (two first lines of (24)) and 3 
(expression (24))), using the cumulants (23). The tails of the distribution are repro- 
duced with increasing accuracy when the number of terms in the analytic expansion 
is increased. Figure 2 displays the same for m = 11. 

It is to be noted that techniques similar to those used to derive (i?") and P{R) in sec- 
tion 1 can be applied to derive distributions for random states drawn from orthogonal 
or symplectic matrix ensembles, since the joint probability distribution for Schmidt 
coefficients is of the same form as the distribution (3). 

The author thanks CalMiP in Toulouse and Idris in Orsay for access to their su- 
percomputers, and Bertrand Gcorgeot for reading the manuscript. This work was 
supported by the Agence Nationale de la Recherche (ANR project INFOSYSQQ) and 
the European program EC 1ST FP6-015708 EuroSQIP. 

Appendix A. 

The aim of this appendix is to evaluate integrals of the form 7(n) given by Equation 
(5). The Vandcrmondc determinant (6) can be written 

(A.1) 



where the sum runs over all permutations on p elements and e„ is the signature of the 
permutation a. For any function ip symmetric under permutations of the Xi we have 

(A.2) 



L 



dx\ . . . dXpV ip{yi) = 



= eaea' / dxi... (ixpx°(i)a;°,(i) . . . a;^^p^)X^7(p)¥'(x) 

= p\ I dxi . . .dxpx\x\ . . .x^~^V{^)^p{J^), 
Jo 
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Figure 1. Probability density function P(Q) of Meyer- Wallach entanglement 
Q for random vectors of size 2"* for m = 10. Top: P{Q) in logarithmic 
scale. Bottom: differences between Pb{Q) (analytical expansion at order s) and 
the numerical curve. From bottom to top on the left axis of the top figure: 
truncation of the expansion at order (gaussian, red, dotted); order 1 (green, 
short dashed); order 2 (blue, dot-dashed); order 3 (purple, long dashed); numerical 
curve from column vectors of 1000 random unitary matrices obtained by Hurwitz 
parametrization (black, solid). 



with X = (xi, . . . , Xp). The integral (5) becomes /(n) = J(T(n)) where the sum 
runs over all permutations of the n,j , with 



J(n)=p! [ da;i...dxpF(x)[|x[+"'+*-i5 ( ) 



Using the fact that 



a\b\ 



dxx^il - xf - -, 

a recurrence on the number of integrals shows that 
nl / p 







dx ^ • • • dxpX-^ 



aila2l ... ftp! 

(ELi«^+p-i)!' 



The Vandermonde determinant (A.l) can be written as 



^(x) ^^e^xl 



(i)-i 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



Inserting this expression in the integral (A.3) leads to a sum of integrals of the form 
(A.5), which can be cast under 



[p^ +rp + }_^-ni- ly. 



(A.7) 




Figure 2. Same as Figure 1 for m = 11. Numerical curve averaged over 100 
random unitary matrices. 



with A(n) a determinant defined by 
A(n): 



1 r + ni + 1 
1 r + ns + 2 



(r + ni + l)(r + ni + 2) 
(r + n2 + 2)(r + n2 + 3) 



1 r + Up+p [r + Up + p)[r + Up + p + I) 
The determinant can be evaluated by recurrence. It finally yields 

nLi(^ + + i 



J(n) = 







which in turn gives /(n) as a sum over permutations of the rii of J(n) 



(A.9) 
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